{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GP Regression with a Spectral Mixture Kernel\n",
    "\n",
    "## Introduction\n",
    "\n",
    "This example shows how to use a `SpectralMixtureKernel` module on an `ExactGP` model. This module is designed for\n",
    "- When you want to use exact inference (e.g. for regression)\n",
    "- When you want to use a more sophisticated kernel than RBF\n",
    "\n",
    "Function to be modeled is  $sin(2\\pi x)$\n",
    "\n",
    "The Spectral Mixture (SM) kernel was invented and discussed in this paper:\n",
    "https://arxiv.org/pdf/1302.4245.pdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import torch\n",
    "import gpytorch\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "%matplotlib inline\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set up training data\n",
    "\n",
    "In the next cell, we set up the training data for this example. We'll be using 15 regularly spaced points on [0,1] which we evaluate the function on and add Gaussian noise to get the training labels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_x = torch.linspace(0, 1, 15)\n",
    "train_y = torch.sin(train_x * (2 * math.pi))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set up the model\n",
    "\n",
    "The model should be very similar to the `ExactGP` model in the [simple regression example](./Simple_GP_Regression.ipynb).\n",
    "\n",
    "The only difference is here, we're using a more complex kernel (the `SpectralMixtureKernel`). This kernel requires careful initialization to work properly. To that end, in the model `__init__` function, we call\n",
    "\n",
    "```\n",
    "self.covar_module = gpytorch.kernels.SpectralMixtureKernel(n_mixtures=4)\n",
    "self.covar_module.initialize_from_data(train_x, train_y)\n",
    "```\n",
    "\n",
    "This ensures that, when we perform optimization to learn kernel hyperparameters, we will be starting from a reasonable initialization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "class SpectralMixtureGPModel(gpytorch.models.ExactGP):\n",
    "    def __init__(self, train_x, train_y, likelihood):\n",
    "        super(SpectralMixtureGPModel, self).__init__(train_x, train_y, likelihood)\n",
    "        self.mean_module = gpytorch.means.ConstantMean()\n",
    "        self.covar_module = gpytorch.kernels.SpectralMixtureKernel(num_mixtures=4)\n",
    "        self.covar_module.initialize_from_data(train_x, train_y)\n",
    "\n",
    "    def forward(self,x):\n",
    "        mean_x = self.mean_module(x)\n",
    "        covar_x = self.covar_module(x)\n",
    "        return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)\n",
    "\n",
    "    \n",
    "likelihood = gpytorch.likelihoods.GaussianLikelihood()\n",
    "model = SpectralMixtureGPModel(train_x, train_y, likelihood)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Training the model\n",
    "\n",
    "In the next cell, we handle using Type-II MLE to train the hyperparameters of the Gaussian process.\n",
    "The spectral mixture kernel's hyperparameters start from what was specified in `initialize_from_data`.\n",
    "\n",
    "See the [simple regression example](./Simple_GP_Regression.ipynb) for more info on this step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iter 1/100 - Loss: 1.339\n",
      "Iter 2/100 - Loss: 1.303\n",
      "Iter 3/100 - Loss: 1.248\n",
      "Iter 4/100 - Loss: 1.219\n",
      "Iter 5/100 - Loss: 1.170\n",
      "Iter 6/100 - Loss: 1.109\n",
      "Iter 7/100 - Loss: 1.088\n",
      "Iter 8/100 - Loss: 1.059\n",
      "Iter 9/100 - Loss: 1.008\n",
      "Iter 10/100 - Loss: 0.953\n",
      "Iter 11/100 - Loss: 0.906\n",
      "Iter 12/100 - Loss: 0.887\n",
      "Iter 13/100 - Loss: 0.811\n",
      "Iter 14/100 - Loss: 0.762\n",
      "Iter 15/100 - Loss: 0.715\n",
      "Iter 16/100 - Loss: 0.683\n",
      "Iter 17/100 - Loss: 0.639\n",
      "Iter 18/100 - Loss: 0.585\n",
      "Iter 19/100 - Loss: 0.542\n",
      "Iter 20/100 - Loss: 0.498\n",
      "Iter 21/100 - Loss: 0.467\n",
      "Iter 22/100 - Loss: 0.402\n",
      "Iter 23/100 - Loss: 0.349\n",
      "Iter 24/100 - Loss: 0.313\n",
      "Iter 25/100 - Loss: 0.251\n",
      "Iter 26/100 - Loss: 0.214\n",
      "Iter 27/100 - Loss: 0.172\n",
      "Iter 28/100 - Loss: 0.148\n",
      "Iter 29/100 - Loss: 0.122\n",
      "Iter 30/100 - Loss: 0.036\n",
      "Iter 31/100 - Loss: -0.020\n",
      "Iter 32/100 - Loss: -0.073\n",
      "Iter 33/100 - Loss: -0.102\n",
      "Iter 34/100 - Loss: -0.163\n",
      "Iter 35/100 - Loss: -0.146\n",
      "Iter 36/100 - Loss: -0.174\n",
      "Iter 37/100 - Loss: -0.216\n",
      "Iter 38/100 - Loss: -0.289\n",
      "Iter 39/100 - Loss: -0.393\n",
      "Iter 40/100 - Loss: -0.430\n",
      "Iter 41/100 - Loss: -0.331\n",
      "Iter 42/100 - Loss: -0.388\n",
      "Iter 43/100 - Loss: -0.504\n",
      "Iter 44/100 - Loss: -0.629\n",
      "Iter 45/100 - Loss: -0.570\n",
      "Iter 46/100 - Loss: -0.578\n",
      "Iter 47/100 - Loss: -0.728\n",
      "Iter 48/100 - Loss: -0.787\n",
      "Iter 49/100 - Loss: -0.186\n",
      "Iter 50/100 - Loss: -0.532\n",
      "Iter 51/100 - Loss: -0.850\n",
      "Iter 52/100 - Loss: -0.914\n",
      "Iter 53/100 - Loss: -0.879\n",
      "Iter 54/100 - Loss: -0.815\n",
      "Iter 55/100 - Loss: -0.804\n",
      "Iter 56/100 - Loss: -0.808\n",
      "Iter 57/100 - Loss: -0.850\n",
      "Iter 58/100 - Loss: -0.939\n",
      "Iter 59/100 - Loss: -1.047\n",
      "Iter 60/100 - Loss: -1.128\n",
      "Iter 61/100 - Loss: -1.181\n",
      "Iter 62/100 - Loss: -1.210\n",
      "Iter 63/100 - Loss: -1.181\n",
      "Iter 64/100 - Loss: -1.044\n",
      "Iter 65/100 - Loss: -0.988\n",
      "Iter 66/100 - Loss: -1.113\n",
      "Iter 67/100 - Loss: -1.085\n",
      "Iter 68/100 - Loss: -1.284\n",
      "Iter 69/100 - Loss: -1.252\n",
      "Iter 70/100 - Loss: -1.305\n",
      "Iter 71/100 - Loss: -1.318\n",
      "Iter 72/100 - Loss: -1.300\n",
      "Iter 73/100 - Loss: -1.312\n",
      "Iter 74/100 - Loss: -1.365\n",
      "Iter 75/100 - Loss: -1.418\n",
      "Iter 76/100 - Loss: -1.484\n",
      "Iter 77/100 - Loss: -1.564\n",
      "Iter 78/100 - Loss: -1.599\n",
      "Iter 79/100 - Loss: -1.678\n",
      "Iter 80/100 - Loss: -1.731\n",
      "Iter 81/100 - Loss: -1.793\n",
      "Iter 82/100 - Loss: -1.790\n",
      "Iter 83/100 - Loss: -1.801\n",
      "Iter 84/100 - Loss: -1.832\n",
      "Iter 85/100 - Loss: -1.916\n",
      "Iter 86/100 - Loss: -1.974\n",
      "Iter 87/100 - Loss: -2.030\n",
      "Iter 88/100 - Loss: -2.108\n",
      "Iter 89/100 - Loss: -2.166\n",
      "Iter 90/100 - Loss: -2.152\n",
      "Iter 91/100 - Loss: -2.119\n",
      "Iter 92/100 - Loss: -2.088\n",
      "Iter 93/100 - Loss: -2.101\n",
      "Iter 94/100 - Loss: -2.174\n",
      "Iter 95/100 - Loss: -2.247\n",
      "Iter 96/100 - Loss: -2.223\n",
      "Iter 97/100 - Loss: -1.789\n",
      "Iter 98/100 - Loss: -2.284\n",
      "Iter 99/100 - Loss: -2.083\n",
      "Iter 100/100 - Loss: -2.164\n"
     ]
    }
   ],
   "source": [
    "# Find optimal model hyperparameters\n",
    "model.train()\n",
    "likelihood.train()\n",
    "\n",
    "# Use the adam optimizer\n",
    "optimizer = torch.optim.Adam(model.parameters(), lr=0.1)\n",
    "\n",
    "# \"Loss\" for GPs - the marginal log likelihood\n",
    "mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)\n",
    "\n",
    "training_iter = 100\n",
    "for i in range(training_iter):\n",
    "    optimizer.zero_grad()\n",
    "    output = model(train_x)\n",
    "    loss = -mll(output, train_y)\n",
    "    loss.backward()\n",
    "    print('Iter %d/%d - Loss: %.3f' % (i + 1, training_iter, loss.item()))\n",
    "    optimizer.step()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Make predictions, and plot the model fit\n",
    "\n",
    "Now that we've learned good hyperparameters, it's time to use our model to make predictions. The spectral mixture kernel is especially good at extrapolation. To that end, we'll see how well the model extrapolates past the interval `[0, 1]`.\n",
    "\n",
    "In the next cell, we plot the mean and confidence region of the Gaussian process model. The `confidence_region` method is a helper method that returns 2 standard deviations above and below the mean."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQAAAADFCAYAAAC7ICzVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJztvXl4G/d17/0d7AsJgiAJ7rtI7SshW3Zsy6bkeG0cy7Tjus1W2Ukb3yRvqip1b3QTt9d25KZ6b5y2N40iN3mvGyeOFbnx9VZbkm3J2ixQ1GKJlESBpLiBIAmCJEBgsMy8fwwGxL7OAKA4n+fRIxIcDA+G8ztzzvmdhaBpGgICAosTUa4FEBAQyB2CAhAQWMQICkBAYBEjKAABgUWMoAAEBBYxggIQEFjECApAQGARI8n0BAaDYav/y7uNRuPfZno+AQGB7JGRBWAwGDaAWfgHAWwwGAxN3IglICCQDQguMgENBoMWwItGo/GbmYskICCQLTJ2AfwYANjiHfDMM88IOccCAjli9+7dRLTXOVEARqPxoMFgeNRgMHQYjcb9sY77+7//+4Tnslgs0Ov1XIjFG/kuY77LB+S/jPkuH5C8jD/60Y9i/izTGMCLBoPhG/5vbQB0mZxPQEAgu2S6DfgLACb/ToDWaDTu5UAmAQGBLJGRC2A0Gk0ATP5vD6byXq/XC7PZDJIkERyIpCgKMzMzmYjFO/ksI0Ewrp5Op4NEwlWIR+BGJWd3iNlshlqtRnV1deCmBQCPxwOpVJorsZIin2WkaRqTk5Mwm82oqanJtTgCeU7OMgFJkoRWqw1Z/AKZQxAENBoNSJLMtSgCC4CcKQCapoXFzxMEQUDo9CSQDAuiFmB0dBRbt26F2WzOtSgCAjcUC0IB/PjHP8bx48fxwgsvpH2Orq4uHDhwAIcOHcK+fftgMjGxywMHDuAHP/gBV6LGxGaz4f77748q1/Lly3Ho0CEcOnQIe/bsgc0WN6dKQIAz8loBaLVaKBQK7N27FxRFYe/evVAoFNBqtSmdx2azYd++fdi2bRu2bNmCJ598Ert27QIAtLe38yF6BFqtFo2NjRGvr1+/Ho2NjdiyZQu2bNmCHTt24Iknnoh6DpvNhj179vAtqsAiIq/3ibq7u/HMM8/gzTffhNPphFKpxEMPPYTdu3endJ79+/dj/fr1Ia8VFxejq6sLjY2N6OrqQldXFw4fPozt27ejs7MTOp0Ohw8fxsMPP4zDhw9Dp9Nh3bp1OHv2LH7/+9+jubkZS5cuxXvvvYdXX30VTz/9NHbs2AEAIcfrdDrs378fjY2N6OvrS0perVYLm80Gq9WKw4cPY3p6Gg8//DD6+vrQ2dmJrq4uFBUVhfysqUmowxJInbxWAJWVlYGItkKhAEmS0Gg0qKio4PT3rF+/PqAgXn75ZQCA1WpFR0cHdu3ahe3bt0On0+Hll1/Gzp078cEHH+D5558HALz33nsAgG3btqGpqQlPPPFEyPE2mw07duxAU1MTDhw4kLRMVqsVTU1N0Ol06OzsxBtvvIHt27fj8OHDAVmDf8YqH4GFgWWWxMDkHFQyMVQyMQoVEpQWyLMuR14rAIDJd37qqaewfft2vPzyy2kFAjs6OvCtb30LTz75ZOC1vr4+rF+/Pqq//fDDDwMAdu3aBbfbjba2Nmi12sDxRUVFgWPb29uxZ88ebN++PfBa8PFPP/00dDomQ3pqaiopeW02G5qamnDo0CGcPXsW69atC/m5yWRCX19f1J8JLAxO90/BNOEIfC8iCHztljqo5dldknmvAF577bXA1y+99FJa59Bqtdi5cyf27dsXMMX/5V/+JeQY1gXYsWMH9uzZg/b2dnR0dGDdunUBE16n08FqtaK/vx82mw1arRYdHR34wQ9+EIhLPPfccyHH79ixA/v370dbWxv6+vrQ1dUV4o6wi/nQoUMAgLNnzwZk6+vrQ2NjI6anpwPHsa+H/yz8vAL5i93lRd/kXMhrFE2jZ8yOtrrU4luZwkk/gGR45pln6OBqwN7eXixZsiTiuHzOsmPJdxk9Hg8GBgaiXt98Id+r7fiU74TJCuNApDWoVUrx5U11SZ8nlWrAWOXAeb0LICCQLD5qYSQ++Sgal0Zno/7M5vRg2ObMqjyCAhBY8FwancEJkzXXYiTFtXEH5tzemD+PpRz4QlAAAgsa25wHR65O4uzQNMZn87/+4cJI/CrSXosDpMeXJWkEBSCwgKEoGu9fssDjo0DTND68PJHXNRCTdjdGEpj4XorCZYs9SxIJCkBgAXOqfwpjs67A92OzLlwYzs8+DUDipz/LpZHsuQGCAhBYkAzbnOgciMzhOGGywu6K7WPnCreXwmVzcgt73E7CkiV3ZlEoAJPJhFtuuQVdXV1xXxNYGNA0jQ+6x0Ej0tx3+ygc6Z3IgVTxMU044PZRSR+frWBgzhWAQiEP+VdYWBDxWrL/YtHU1BRIBGKZnp5GY2OjkDyzALHMujHr8sT8+bVxB8ZmXDF/ngtGplOT56rFnpV4RkaZgP6BIOxosI35PBosOH03HJPJFFLAAyCiCGfPnj14/vnnsX//fjz77LNZklogGsnslQ9OuVCuUWRBmuQYsaWmAFweHyYdbt7rAzK1AB4DoGNnAQS1CE8al4sM+Tc7a494Ldl/idi2bRsOHDgQkTa7a9cuNDY2orGxES+//DKampoCacBvvPEGtmzZElILIJBbhpNYTCPT2U2oicec24epOXfK70tVaaRDpl2Bg9uAN4FpE563bNmyBU888QR27twZ8bPgAp5oRTjFxcXZFFUgBhRFJ9xKA4BRmwsURUMkyn3budEUzX+WkWkX1tTEtly5gJNiIP9QUKu/TXhMLBZL4GuKouDxRPpxPh/3SRBnz57FL3/5S6xZswbr1q2DWq3G6dOnYTKZcPr0aTz77LN47bXX0NjYiOLiYvT29qKurg5WqxW9vb145ZVXYDKZMD4+DpPJhM7OTrS1tXEuJ1f4fD5QFBVyvfONdC2pcbsHUzOJA2QOAD0DIyhVp1ezwaWld2lgFg6HI/GBYVwZcmFDWWwFxoWMXFUDdiQzGDS4cGFmZiZmQQ3XhTYbN27E7373OwDA97///cDrJ0+eDHzd2toacjzLo48+CgD48pe/DAB49913874YCABEIlFeF9sASEu+604b1OrknqikRA29Pv3qOq6u39x1N9Tq9N4rKyiGVhX7XstUxox3AfzzAP/R//XWRMcL5DdOtw9nB6dzLUZMUimWSSZWwDduL4Vxe+r+PwvfsYxMZwNuBfCiwWDoNBgMnRzJJJAjHKQXb5wdxXHTJJzu7OWjJwtF0Sn50yM2V85Tg0enM5OB70BgpkHAgwCaOZJFIIfYXV7857nRQLT64ugsDPXZbU6RiLFZEp4UkmlIb3a20uKR6v5/OHxbMTlPBBLIPdNOD/7QNRKyVfXZ8DSoPKuxT6dWPtduQKYKYMbl4TW1WVAAAjjaO4mZsMy6WdIb0rMuHxhKYzHnUgH4KJqTjMThDJVIPAQFsMihKBrDU9GfrOfyqLLOR9EYTWMxJ5MzwBdjMyQnnYr4/Ax50RT0nz+8Fvja56MgFqeul759V/xQRFdXFzo7OwONNDs7OwOtvROxb98+tLW1Yf/+/SHvsdlseOKJJ/DOO++kLG++MDZLxixSGbE5MT5Loqwwdz40y9gMCS+VvP/P4vT4YHW4oVPLeJAqPlxF8DN1I+KRFwqAb2w2G37yk5/g1VdfDbzW2ZncpoXNZkNfXx+efPLJiMk+sab9LCQGYzz9Wc4NT2PrstznE2TSK2/Y5sqRAuBm4U45PHC6fVDKxJycL5hF4QLs378/YgTYzp07YTKZAvMC2f/vv/9+dHV1BeYFdnZ2Btp2syO72FFjhw4dCrTqNplM2LdvHw4cOACTyRT1XACwZ8+ewO8Kf08uSKQArozZ82JLcCgDBZCLugCapmGe5qamnwbNmxWwKBRANLRaLXbt2oX29nZs2bIF+/fvj1r009bWFpjdxz7tf/jDHwbex74WXlAU7Vz79u1De3s72tvbMT09HfGebOPxUTAnuLF8FI1LSTay4AsfldliGp7KfiBwwu4G6eVOcQoKIAM6Ojpw+PDhkNfYQRzhJFv0E23aT1tbG9avXx+IE4Sf6/Dhw2hsbIRWq8W2bduiviebjNhcoJJIUsl2q+pwpubccf3/uVkRpiyxvVmH2wvbXOz+AXyQbgFQLPgKBC6KGED4ZKDp6Wm0t7ejsbExMMWno6MDXV1d6OvrC/j9XV1dsFqtIa+ZTCZ873vfi5j2Ez4RCEDEuXbu3BlyTPh7st2cJJH5z2KZyW233Yk4qbQ9RjVefbEGXg+BHT+/hpLK6At9bJaMm1PPNVy39Bq3u+H2UpBJuH1m54UCCI7g81VoEzwAlEWr1UZM1WUj+sGR/S1btoS85vF4AsVDJ06cCByXzLnCZcjlVN9oCmCoV4HuTwtw2xesUBYwT12nx4dZlxeFitzcLtEUAE0DH71egnd/rQdNMxVzH75eio7vjEY9x/gsiaXlBbzKGcwYx0qTppmcglqditPzLgoXQCASp9uHyaCFZbeJ8fpPK/Gz7zTi/Vf0ePfXoZH/bDWpjMaEPfR3u10EfrO7Gu/8qhw0TeCWB60gCBrGD4owPRFdSY3bsye/20thigeXg484gKAAFilDNido0PD5gI8P6PDi9iX49L+KQfjviDOHi0A652vRc9ljL9wC+NWzdTh3pAhypQ9f/eF1bHvajNW3zcDnFeHjAyVRzzE+m35FXqpYZsmoDUszheu4ApBDBUAQRM4rtW5UaJoGQcTvhMOa/4d+V4q3flkB15wYywyz2PFv19CwYg6kU4yzH893o8mVBWAnvXAGTcqxmqXoPaeGQuXDt/9XH1bdwgzRaH9sEgBw8p1iOKYj98tJrw8zzuwEAsd4ulbmGZLz+oycKQC5XA6bzSYoAY6haRozMzOQy+Nn7w1NOUHTwJlDTMXfl/56GNv/5yD0NW5sup/Z2Tj5zvwuhmXWnZO/1WTY0/9yJ+PHt6x3oLx+/mfVS1xYtnEWHlKEo3/URT1XJnX5qZCMsvR6mDhGKnh8mfUWiEbOgoAVFRUwm82wWq0hNxZFURCJ8tszyWcZ2Sd/fX19zGNmnB5MOz2wDMowOSqDSuPFhvb5JiBrbpvBm78ox9BVJYauKlDT4gLp9cHm9KBYld2MuvAb/nIn01pnqSFyfNaWxyfQc7oQx97UYfMjk1CqQ7cOLbMkmsvSbM2TAokCgGMDMvzz9xpRWOzF7Q9bsfFuG6Ty5LTB6LQL5RruUrNzpgAkEglqamoiXs/3ufFA/stosVggkcT+07Lm/6WThQCA5RvtEAVZzVI5jbat0zj6RglOvlOMju8ykXXLrDvrCmDSMb+YfF6g96xfAbRFKoCGFU40rXbAdEGNE28XB9wClmwEAp1uX9yZBQBw8HdlIJ1ikE4x3vjXSrz/Shlu/YIVtz80v/MSi5FpF9bVctcoND8fYwK8wpbIXjzFKIAVmyIz/Tbdx7gBXR8VweVgbpNcxAEmgoJ3A90qkE4x9LUktGXRa+S3PM5MBTpyoAQeMjQOko1AYCL/f2JEinNHNBCJaWz7b6OoaXHCMSPBB/+hx/95rjbh+bkOBAoKgGfyMcZhmSVht4lxvVsJsYRC64bIun99rRtNqx1wu0Q482FR4H3ZxOujMBUUuAuY/1Ge/iwt6x2obHTBMS1B38XQPfM5txcOkt+5gYmu0Uf7S0FTBNrabbjlgSl856U+/OWL/ZArfeg9p8b4cHwLa47jrEZBAfAETdO4PGbHK6cGs56GGg+3l4JtzoPu0wWgaQJL1s1BoYpudt7CBgPfLQZNM8k02ewSNOnwhChQNgDYGkcBEASjBACg77PIpBm+A4Hx/P/pCQmMHxSBIGjc6XdPCAJoXjOHVbcyVljXR5qEv4NLK0BQADwwYnPh9c4RvHNuDKePy/Dup/nTWGPS4QYNOuD/r7g5dqHPqltnodZ4MWpSYPCyEh4fPwkusQj2/+02MYZ7lZBIKTStmov7vsaVzM/7LkVRADxbMfEsgCNvlMDnFWH152ahrwlVROvvYoKwXYeLEu4OcJkQxEVb8K0Gg+EDLoRJhtHRUWzduhVmsznk63zh2rgDP/+/43j5pxo8/5VW/PK/12PnX5TCMp29RJR4jM+S8LiJwNM0ngKQyGgY7mYqGc8dZZ5M2XQDgp/WV84w5n/T6jnIFPFXSMMKRgFc71HCF2bx8xkInHV5MeeO7mI4ZsSBbdW7vhQ5vXjJOgcKir2YGJFj8Er8mYZ5ZQH4OwNnjR//+Mc4fvw4XnjhhZCv8wGHA/jy4wr85Kkl+Hh/Kew2CUQiGhMjMvzT/84PN2Dc7sa1c2p4SBGqm50xg2ksS9Yx5jR7U/KV5BKN4ByAy2cSm/8sBVofympIeEgRhq+FLiYLj4HAeMrx2Js6uF0itLbZUbMkcgGLxcD6zX4r4MP4Uf6pOTdnPRqyug2YzKiqWOOOWlpaQJIbwcwj/Rn27iUAEAD+O/bu/TH27t0LuVyOq1evcilySjL++69VOHOsBFKFD6s+Z8WGLeOwmhV442dNeOXfVPjqoyMo0/B/yeONjDKNTOLcJ1UAgCUbrAlHVumqSQD1GOpVYHbGgV64sUKbemuuVGRk6Tdb4fZRoCjgspGxAOpWjMPhSPwErFk6g/GhMlzukqCkZv4zOhzA4IgM8gRVdemM3bo8GH0EGOkU4egfmaf/rQ8NxrzmyzaZcfQ/S9D1kQZ3PWGCOE4DoM/6hlFEZG4JZFUBJLt3Hu248+d7sHZtEVyuyHp9qdSMRx6xY/fu3Vnbn4/2e/b/gbmpOr49irsecKJCUwidSoZP/kBifFiBX7+uwU/+LjsVadHk81E0SMKOK53MNVx3Bwl10MwqAgS0KikIAFZ/i3C1Gigud2NqTAaHVQe31o2S0jKIORi6Ge9vNevyQqqYhRTAcK8Cjmkpiko8aFgqBkEkTuZpWetG1yFg+IoWanWY1aAsgr5YmZF80SBHvFFHgF35VAOXXYK6ZXNYYaBjyt+yBiirJjE+LMfoVT2WtsVWzqRYDa1GkfvRYNni9Olq/+K/ConkewCeBvBzAIDHsx0ajQYVFRU5k6+7m8Bn56RQqin8ZIcOX9lUh8+v0MPQoMVX/or5Q766tzCnsQCrw43Bq3LMTDKLqbqZeYKsqtJg2/oqfOP2Bvz5zbXYurws5H21rUzi0OBVBXwUjUkH/58h2Fe/7Pf/W9vsSFDiEIANBPZfVEUE1fgIBNI0HTPP4Np5Rv5Vt87GlZ8g5oOBZw7HdwO4igMsCAUwOjqKp5++AAD43OeMOH78CTQ0vIO6up9BXeAFsAkXLkhyGhD8918z/3/hix6UakNtt7/+pgylVW5MjsrwT/+WOwUwbidxyZ/8s/xm5mYUiwjc2qRDtVYZaDZRrlGENNGsaWFutsErzFOT61r3aAT7/1f8Act4T8RwSio9KCz2wDEjidhb5yMQODXnidkCzHSB2Y1oXh1/9wKYVwCfHdfA7YqtLbj6DJwMB2X+M3RwIE9Uvv/9X2N6egMkEhfeeGMb1qxZg56eHly5chZf/xpzzNTU4zkLCHq9wG9/yyz6v3wy8o9WUijFl7/JRNtf3VuAOVfmPnQ6jM+6Azfjso2MWdxQooJcGulsrqgsDHzNWgBDVxkFkI2UWnYHwOMm0H9JBYKg0bI+cQCQhSCAxpWM3P1hCUF8ZATGCgDOWCUYH5ZDpqBQvSRxB6bSKg/qls3B7RLh4snCmMdxlY7BxS7AfqPRWGw0GvdzIVAwWq0WCoUCr7/OVHd5vf8f9HoFtNr5mXV797YBAHp6NoCiCrB3714oFKHH8M3BgyJMWMSoa/Ri06bof5m//isZSirdmBiWY++vc9Nld2yGxNBVJipet5S5GZeWR7/JlpYXQuS3V6v9UeuRa3J4PZEVenzANgEZuy6Hz0ugtNoNVWFqirOBzQcIUwC2OQ/cXm6VsDmGVcQq3IYVcxAnGXHbcFdyuwFckNcuQHd3Nx555CsAvgoAkMt/hccffxw9PT2BYy5f/iP0+s8AqAF8BUqlMuIYvmHN/699hYrp45VpZHjiG4wV8G//kv1BGzRN41IPQDrF0JZ5UFjsg1wiRr0uejBMJROjsZS5eZVqCmU1JHxeEcz9CiaZiMcUZ7eXwoyT2Z4c8W/jVTWl7vMGEoLCMgJpcB/HiDXFl1UATUmY/yxrbmcSx66eVcPL8+5xXiuAyspKWCxbABSBIE7C4zkdEeyrrKzE6tWfAAAI4ltwucisBgStVuDdtyUgCBpf+XL8RfE3fyWDXOlDf68UQ0NZES/A1JwHfZcYxVPjN+mby9SQxJnCtKxi3jqoaZkPBHp8FGZ4HFhp9WcrAsBoHyNzOgqgsskFudKHyVEZZqyhj18zhx2OSP/0oWjMK4Dk4xeFxT7oa0l43SIM9yberciEvFYANA2cP/85AMCzz+rx1FNPYWxsLOI4leoDqFQ20PRyPPjgT6Iewxevvy6Cx0Pg9s1eRKluDqFcK8NqA2Mqvpe13EmGCbsbg34fnvXpEzXJbNCpoJJJ/O9hFsyQPxAYr1NvpkwELaYRE2MBVDamHncQi4G6ZWwcIHQhjXI0tAMARmeitwCz28QYu66AVE4FrnmyNMZwX7gmrxXA6dMEpqebodPR+O53q/HSSy/htddeizju9ddfxXe/y9zMcvl3ox7DF7/6NXMJt389ueMfvJc5/t33+ZIoOuN2MhDFr21xoVAuQbU2fsqpSERgWQVzXectAOYcsZ54XMDGGGh6XgFUNaf3xI61kLjscRhrS87kdz3qlzkhSbHRdeOq6O4L1+S1AvjFL5jo9Ne/7oMi/r2K7dt9EItp/Od/ijAavTM05/T0EDjbJUZBIYUvfCG5oNJD9zOX/PhRScotoTJhdMqNkWt+F6DFidbygoR9AwFgud8NqG52QSSiMdYvh9tFhDyluYY9t80ihcshhlrjhUYX3eWQikWo0iohl0RPmwvkA4QVBtlJL2Y5cmNiFeekY/6zBBTARRXSmImaNHmrALxe4M03GfH+4i8SR81raoB776Xg9RJ4++3sfKz332d+z5/8CQVlkq7asmU0dGU+TE2I0d2deTZdspw7D/i8IpRVk1AWUCH+fTx0ahkqNArIFDTK60lQFIGRPgWvOwHsuYdNfv+/2RUSXCUIAu1Ly/CnG2vwjdsa8Mj6KiyvjP556pY5IRLRGL6mgGsu9L7gIpnGR9ExrQnThfkCplQp1ntQVOqB0y6G5Tp/QeO8VQDvvz+J2VkC9fUeNMef/B3g7rsZVfnssx9nJSHoyBHmrmy/K/lHOUEAt9zGPHnefi87JsCsy4trQQHA0gJ5StNyW/ShbsDQFSVsTg+8McaKZyorm1AzGsP/r9EqsLJKg9ICOUT+lOSVMRSATEGjqtkFmiIweDksDsCBG2CZJeGLsik/NyuCuV8OiZQKxCFSgSDmrQATj25A3iqA3btPAgCUypNJv2fzZuYPMTGxCs8/z29CkM8HHD3KXL7bb09tIdz/eeb//8pSHSXj/zOLqbbVlfKEHLYJJRsIHLyiAE3TsPLQGyB4CEgs/7+xNDKXnrVUosEG4MIrA8c4CATGMv/7LqpA0wTqljohlaWn6LMRCMw7BdDS0gKFQoFPP2X+yD09/zupxB5mEq8cgBlABX75yyO8JgRduEBgeppAXR2FhobU3nvP3cz/xpNSeLJQJTw+GxQAbHWiNolCmGBKC2QgCGLeAuAxEDgZZQegqjF0kTVFUQBAaPZiMGwiU7gCGLeT8GRoxcQMAGZg/rM0BQUC+YoX5Z0C+OSTT9DR8QSA2wAACsWppBJ7uru78aUvfQli8VEAgFR6D68JQezTn7U6UqGmBqhr9MA1J8Knn3ItWSSD4x6MXZdDJKJR3+JGSQrmP8AE2krUMlQ2kBBLKIwPyeByiHjZCmTP6XKIYDXLIJZQ0NfOP6lLC+QxZxS26AsgjZLXwBY9DfeGKgCKpjNqcELTdBwFkHoCUDj6OhLKAh+mJ6SYsvAz2DTvFEB5eTlIcg0ANQiiB273QFKJPZWVldBoNKAoZgy4x3MbrwlBrP9/xx3pPUHuuJPxc9/6L/7jAGfOADRFoLyBRK1eFvCbU0FfKIdERqOykQRNExi+puClKpANALIJQBX1ZEgKbVNpbHNYJhEF4hXBVNSTEIlpTAzLQsadAZnlA0zNeeDyRAaonQ6mEYlITKN+efoKQCSa727E13Zg3ikAALh6lcmo2bZNFzP5JxoWiwWPPcaUsioU98BsTtyAJB0oKn3/n+WBe5kb8dAhfv8EdtKLq5eYJ35tqzOmn5wIfSEbB/DnA1xRcr4TENwFeKTPHwBsCl2g0fz/YKIFAyUyGhX1LtA0EXArWDLJB4jl/w9cUoGmCNS2OhO2L0tE8HYgH+TFePBwKir+FD09wCOPFGPbtpeSft9rr70GmgaOHKExOlqIH/7wNYCHIY3d3VLYbARq0/D/WdrvBAgRjYtnpbDb3SjgqU+IJcj/r2nJXAGwpcHD1xRwuCfhdPuglMVpXZMC1rn5LsAB/z8oBbhQIQ3IEYuKIgV0KlmgoQlL9RIXRkxKDPcqAlWCQGZbgbHeO9DDXG+2GCkT+E4IyjsLgCSBkyfTN68JAti8mXnfxx/z8/FOnGBuws13pK9cioqAFas98PkIfHyEv3yAsRkSQ/4dgLpWMm0FUFogg1hEoKKBuenN/cw14NINCLYoRqMogMY45n8wK6oirYD5QGBoANTp8aXdtj1WAdBQYMcl9e2/cGqWOCGVU7AMyqMOPc2UvFMAXV1yuFwEVq2iUFqa3jnmFQA/C+vkSebmT9f/Z2lvZ97/Do9xANOwBxMjzH70shXpP63FIgKlajkq6hmTfHxIDp+XWwXAZgD6fMCoX8FUBu0ANJUkpwCWBZUyswQUQG+kAkynMMhOejETZQQYTSNox4WDnn3S+dJtPtyAvFMAx4/7n66b019c7HuPHBVxnkZJUcCpU9wogAfuYS7/hx/y92c4c4ZZCFVNLtToMsso02vkkClo6Crc8HkJjA/JObbG9n4CAAAf8klEQVQAGOUyMSyD1y2CVj/fA0AuEaNam9z2pVImjqhzqGp0gRDRGBuQw+MOVQ6xavnjEcv8n7JI4ZiRQK3xoljPzR4vn/kAeacAAuZ1GttrLI2NQG0tDdsUgfPnubUCmP1/cVr7/+HccgsNqYyG6YoUVisn4oUw4/QEMgBrW12oKErP/GcpLwx9Kpv75ZwGAtktwHnzf35h1pcoU9q9KA9zdWQKGvoaJpWZdV9Y0okDDFij+/eDQeZ/sv0LE8FnHCCvFIDTCZw5IwdB0GlH1wF+4wDz0f/MzXa5HFi5lnlKHD/OvbsyNjufAVjT6kRlmv4/CxuAq2hgFuZov5yz5iB20gunf0stWgJQrOSfWFRGUXax3IBJhzulDkEeH4VrlugFPoGAKwfmP0v98jmmnqFXEbdPYDrklQJ4++0peDwEVqxwoziy+3dKsApg9+5TnNYFvP8+c/OvXTvFyflu+xyzeA4f4T4OYJklA/X7S1a4oVNnlkxSrJJCKhYF4gDmAe6agwRbEvNbgMwiIggC9brUnn6VGjkIhC4WNqU4PCOQpmMX9ETj2rgD7hgZhEOBkuvMA4AsciVTz0BRBK73cNsghJOmoP7xYN/P9Fx79pxhhBIdyfRUAQUwNbUKzz23O+PzAYz//9FHzNdnzvy/nJyz/Q5/eTAPFkCPyYvpSSkUKh9WriCSKv+Nh0hEoLRAhsqwnQAuMgKDy4tHTaFdgIpV0kDH4mSRS8UoVoUqvJqABRC5iC5bkm842mOOPk6NohDouVjDwQ5AMGwcwMRxHCAjBcB2AvaPB7MZDIat6ZyHbf7Z1cU0Qbxw4aWM8vi1Wi1aWxUArgEowr59pzOuC9BqtVCpboLbrQYwgN/97sec1Bps2sQoqksXJHBzmFdD0zS6zjAR/+oWF6q13JSUlmsUKK12QyyhYTUzmXVc1ASwFoDdJsaMVQq50gddBeMelRWkJ3t4zINVKKN98oiZgVfG7JhLYtzWrMuLoano1sL4kAykU4yiUg80Om4bvzYEzTngkkwTgTYCYNvvmABsABCzxi3WaLCjR4/i2Wd/grffvgmAD3L5adx33xexa9eupMaJRTvfc889h7feOgKfrxkSyefx4IO1aZ+PPefXv34On30GAB9DoVDg3nvvzeicLLUNegz2y3Ho0BTa2jKPHNtsNticXlw+7/fZG2cg8zpgsWRuqks8TrhIB0qrnRgbUKG/h0ZN0TjqVanJHT56q888CcecB33dzB6+vm4OTifjZ4s9orSusdznDB3DRQDFFS5MmRXov0yhoiH0Kf3JpQGsqyqIKh/L2RE77I7o1kLvBdZ1sSccuValkWOG9MJOJqco9A0kgFr0dysxM+OARELAZuPg75nh+8MffyXxDo41xkiv16OysggEsQJi8QZ4PBPQ6/VYtWpVWkLp9Xro9Xp/XcDX4fXeDr1+Mu3zseecnWUul0TyCdxud0YyBnPLrcBgP3DmvBb33ceNK2CllBjz+9LNK31Y0VAVtf9/qkgL3DCavahq9GBsALCNaeEk7GmNqGLf46NoeEV2qNUyTJmZW6q62RMYW9ZaWw59ijEAAJCo3egKU3q1LSSmzApMjRSjeWWoATw8J8bW0rLAbkO0zzRmcoWMUwtmfICRvXGFO+YxAJNT8cWNtZCKRXjn4hhGbIndBbUaKK0iMTEix8xYCWpbSWi1BTkfDWYDoMvwHACA8XELvvGNLXjrrb9IKf8/FhaLBX/2Z8wQTKn0LoyORo5kTgWfDxgeXgIA2LfvzziRkeWuO5gb7ugxTk4HgMkAZCPSq9d5OVn8AKBVSiGXiFERtBWYaXutSYcbFM12AY5sAlKapgtQrJJCEfa5Y5UGA8xORO947Cf36LQLNmdsS2d+xyV+QHFtTRE0SimUMjG+uLYSK6s0cY9nidXmPBMytQBOY94KaAKQdq9btpGnxWLBnXfemaFY8+c7dYrC1asK/M3f/BaZ1AWcO0fA6y1AfT2F9vYmPP548jUKibidqXxGl1EKmvZysn98occHp12MgmIvVrRwV/JBEAT0hbL5nYB+5qYfnXahUJFeQUPwPjxbBcgGGtUySdrZiwRBoEKjQP/k/KKeLw2OHk0/NzSN1hgNU7pjBP8AwOuZ376MtwOgkIphqJs3nMUipr2ZVinFsWuTsT8MgIaVTpz+oBh9F1W4Yxs3u1AZWQD+aUBNbPDPHwzMK9iEokzzAdj3Z5KgFIvmZhpFxT5MW8UwmTI/H0XTOOcPANa2OFFZxG1POX2hHJUN7FYgc+5M2msN+3PqfT5gbICNWzDnLy1MrXdBOOGfvTpoKzBalqh5xhV19qHXR6E3xt4/AIwNKOD1iFDq77kYi5saiqNaY+tri6BMYKUFLIBL3DUI4WI02D8ajcaDRqNxLxcCcQ1XCUFsXUEmKcqxIAhgw0bGhD74Uebntzm96O/xP42Wpl8BGIuyQjm0eg/kSh/sNgnsNjHMHFTVTQzL4PWIUKx3Q6lmrkO6OwAs4RmBBVoftGUeuF0iTIxEVy7nh6cjXrs2MRdz+CcQlAHYEvs6aJVSrIph7hMEgeay+MlOpdVuFGi9sE9JMDHCTYOQvEoE4gM2X//ECSLtbTavFzh2jLUA+OnRzCYEfZh5CgQmHN6A/9+0jIzYD88UnUoGgpjPCDT3yzFud6fVXmtqzo05N6P8ovv/mVkA5YXyiPwHdkgne43CuWqxw+mh4PL40GOexdsXzDjcMx7394QPXYnGrc0lEMdJZ16SQAEQROx5h+lywyuA8nJg+XIKc3METp9Oz7k+c4bA7CyB5mYq4fSfdGn3BwI7P83cXx+b9gSCXOvbqIwTgMLRqqQQEURQSrDCn02XelFNcEmtOUoFYKYKQCYRoTSsBRrbpXcwRladj6Lx9iUr9h0bwAfdFpgmHPAmqCoL7rkQjQqNIuETvlqrjAhahhOYc8BRIPCGVwBA5nEAPv1/lg0bAKmMxmBf5oVB5y+J4SFF0FW40VLLfS85sYiAViVFRX1oRmA6ZbUjIQHAUAtAKhZBq8xc/vC6gHq/AhiIk1Y7Q3qTrnFwuwiM9ctBiOjALkM4iToZAUymZaKah0a2RdhFblKCF4kCYLT3Sy+dS6su4P33ma2ftWt5KNnzI5cDy1Yzv+fgx+lnkdldXlz6bL59V7gPzBUlallgobIKIJ3+esEWQGAHwG8BlKhlnFgv4RmBNS1OECIaIyYFPGTm5x8xKUBRBMrryJgtwGoSjGFjSeQGVC1xQaagMD4sx/h45st3USgAtrJwenol/uEf/jGl97rdwIkTzGU6fTq196YKmxb8/ofpxxn6Jucwco25iRgFwM9UmRK1LKQoiKJStwDsrvmmGk6HCFMWGSRSCiVVTLCmLEH7r2SpDLsGciWNigYSlI8I5O5nAjtwJJb/LxOLErYyY6kpjj3mDGAHnjJWgNGY+fW54RWAVqtFTY0CwDkACvz7v19MOodfq9VCo2mH1ysH0I3/+I9/gkKhQEtLCy+ybvHHAT49lX7STt+kA8NX/Rl0qzwokPPT9lGnlkFd5ENhMRNRt1mkcMUZkx2NYPPf7H/6l9eTEPs/fqb+P4tGKYVaFnod6v2LaKAnc1+aLdBpWBFdAVRqFUn3MhCLiLidjwEEehqePi0ogISEzwuQSD6f9LyA7u5urFr1Hf93H0KpVOLxxx/HsWMcpuwFwSYEmbrlsEylXhPg8VEwmV0Yu64CIaJhaOMvZsHOFgjuDQCkZgWE+P/9kTsAZRwpAAARlhAbCMy0vJam5zPz2MYd4dQk2cmIJVGwkA0Enj6d+fW54RXA/LyAQwAArzf5eQGVlZWwWtcAAKTSYyBJEhqNJuP861iUlABNS93wukU48G7qCuC61YnBXjkzA6CORJ2ev6GSRUoJJCJR0FYgs4BTaa81HJQDH+7/EwSR0vzCRIS7E/UcKYDxIRkc0xIUFHtRWhXd+km2lRlLnU4FWZQBJ4GfL2MahHz2mQz25KuYo3LDKwCASS/+6lcbQRA0RKJbMTycXBolSQJmcxMA4K23/obT/P9YbG5n9sTf+6/Ug1P9k3Mh/ihf/j/ALlBpIGWXXcDJttcivRSmHPNKLrAD4FcoWqU06pSfdAm3Jkqr3czUnUkpbOPpu0ns079plSNqCrdcIoY+xWxGsYiIu2sgV9K492vjeOmlyYC7lC55OReAa9i6gHPnaHR1SfCtb/0HkqkLOHWKAEXJsHIlhc2bV2DzZib/P9Py33g8cB+NX/0cOHNcCbvLi4IYY7DCoWmaUQBXGOukttWF8kJ+5iKylKhlqG5hB4YyimfK4QHp8SUsPhqbdYP2/w0oKjIHgCv/nyXcAhCJmCzJK50FuN6jhLYsdp5/PFj/vzHGDIAqrSKtnYzmMjUuj8WWqf2xSXxhRQGUyuTGvMdiUVgALFu2MNH1119PTm2yx7Hvywafv1MMucoHy6Acx84l31XGMkvCQXpx7Tzz5Fiz3pNyF51U0allKK8jIZVTsJplsNvEoEEn5QaYZ+ef/jaLFKRTjMJiDwq0zBZouhWAsSiQS6CKCASy+QDpBwIT+f/h3YmTpbZYyXkCVzQWlQL46leZhfy710SI0e8hwMwM8Opvmcvzta9lTwHIZATW3sQ8Bd98O/n3mSbmMNyrwMykFJoSN9rW83/zlKhlEIvns99YKyCZwqCxkB6A7NOfnwBg4JyF4RmBzKJNNw5gG5dgakwGhcoXcF3CSTUAyCKTiHi5BuEsKgXQ0kLjrrt8cDkJ/OY38a2A3/5WBIedwO23U1ixgv8BnsHctYV5Cp76RJ5UmyqA8f8vnWLMwdY2W8YtwJOB3QlgB2CwCqB/Yi5uFp3bS2HSMd8/YD4DkLsU4GiEFxaxAzeGehXwptGIiX36N6yYgyjK7SSXiDP6HNE6G3PNolIAAPDNbzJP8127hjA6Gj0rkKaBf/1X5gZ+7LH4Ndp88NADzNP7apcKl0fit5YCmD51E3YSl04xdeytBhuvAUCWAoUEcok4sJDYirhxO4kec+zwdM/YbKABCACY/QqA3VFQySRQ85C/EJ6MoyqkUFZDwusWBWYRpIKJNf9jjACvTtP/Z6kSFAD3PPggBZVqCk5nA7797f1Rjzl2jMCVKzIAZpw9+6PsCghg3TIpKupJkE4x3j2c2ALom3RgekKC4V4lpHIKS9bMRhTA8IVOLUUtqwAuKwN16sdN1qi99mecHpy4FppSHZECzJPpGy0bry6JuoBY9CUIANYUZ7bFKFgAHKPValFQoMDc3E8BAG+9VReRFajVarF16x/83+3Dyy//nJPuv6kgEhHYeBuzGI5+JEnYbss0MYdLn/qf/hvsqNBKUpqikwklahmK9R6oi7xwzEgwNcYU78y5vfi0P3K79fDl8ZCe+m4XgYkRGUQiGuW1TFyAL+VVqJBEVNvN5wOkFgh0zIgxNqCARErF7AGQbgCQRS2XoIiDYqh4LCoFwGYFKhS/AeAF8EU89NBfhWQFHjlyGQTxKAAfgL2B7L9kMge5ZOtWZpH0nC7A2xfMMSfXXBiewaB1Dt1+/3/FzXaUZTgAJBWYgp15f/r65fmn3rmh6ZDU4M9GZjA4FbqzMXhFCZomUF5PQiKjA+fki/DtwLpASnBqT2u2Gq9uqTMgdzAKqZiTz8G3G7CoFACbFeh290EkeguAFEND94ZkBb77bhVoWgLgLSgU44Hsv2QyB7nkvq1iSOUURvsU6B3w4WCPJSKwdm3cgY+vTsDtInD1LLP9t2zjLEoLspfeoQsEAiObbFA0jaO9TAxl1uXFsd7IeAq7bdm8Zj7WwZcLAAD6sHNXNIRuYyZL32eM3LG3/7jZxuPbDVhUCgBgknieeuop/OxnywAAn312C7ZsuQdmsxlDQ6N4/gWme/D99w/gyJEjWcn+i0ZNiQyt65mb67KxANfGHThhmvedR6ddeP8SoxSunlXD6xahdqkTRToKlRn20UuFiJ2Ay6FP0uvWOZgmHDgUZvqzXDvPmN7Na5jPShAEdBx3MAom3AIQi+eVVypddgL+fwwFwFUQlm8FsCgyAYNhswJpGvjpSxR6r5bh2LHXsGyZHBQthpssRKHGgv37vwmRCHjpJe66/6aCSERg0+0kLp4swOVONW66x4bO6zbo1DLoC+V464I50KXmUsD8n0VZoQyqNLvopoNSJoZKJgksoqFeBXw+hKSovn/JErVdmMfNzLojCBpNqxkLoFgphYTDFOBwovUYXGqww3RBjXNHNFj9ucQZgW4XgeFeBQgRjfrl0ZO1yjkqZdapZVBKxYHBqVzDxWzArQaDIe124LmiuFiL3qtP+r8rh8ulhZtkFtLszA+gUmU38BeNe+5hTP6rXQWBqbCHesbxx3OjcPlvCIoCuv3bfytunk2q8wzX6NRSqIt80FW44SFFGLseevPH6hV4vUcJr0eEikYSqkLmGD7Nf4ApYgqvt19/J9ME9OKJQjgdiZfEQLcKlI9AdbMLClXkZyNAJF3/nwx8WgEZWwBGo/GgwWD422SOTSaHPtZIJq5hx4e99141SJIGQTDDMynKAYXChXvvjT2aLFsyLqn2oLJJiVGTGu/+nyJs/bNhAEDwM2q4V43ZKSk0pSQ0eisKaAI2W+LcAS4Re+bgcMyhqnkWVnMJes+LUKSPL4PL5UK3kTH165bbAqO00h0DlgoyyhUSnJSpgYaVM+i/qIHxkBwbtkzA5YqdzXj+GDMLp6Z1OuoIsCKFBDZrZoNoglFQzojfQxD5MRosJZIto+Wr3Db8d+j1eng8ZigUMrhcLtA0oFAokhr7lQ0Zy2gaj393DD/9f1Q4/mYlbrrbGZIuCwB958oAACs3OVBRUoQVjdWwWCxZkY+l2avAgH0cjSs9+OwYYOnXQq1OXMcw2M1YWMva5kdpLanRQ8+zFdM8I8bsYKgS33i3Hf0XNbh4TI/bv8DIHm28F+kkcPZD5prffO9c1GOWVBRyev1Xyl24ZA21NAiCyIvRYAsaNiB45MgRNDQ0oKGhIaeBv3AIgsCttwCbHpgC5SPwh3+uDBlmMdonx4l3igEw5n9DCbeTY5OFLdxh4wDXY7TbDsZDEhhg/f+gQBrXRUDRiJZjv/q2GUhkFK6dV2PKEvu5aDyohWtOjMaVc4Fx4+Fwaf4DTOBSIuJnqSa0AAwGwzeivGzKxylAqcIGBAGE7PPnKvAXjWqtAvd9zYKLxwsx0K3CqXeLccsDU7jeo8S+/1EHp12MlvV2tKx3oLG0MicylhbIICIYn1gkojHWL4fbRcRskAkAQ1cL4POKUNXkDPj/cokYhUmWP2dCtAWqVFNYuWkW544UoevDItz8YORwEIoCjr3JmP+feyh2ijhXAUAWsYhAuUYe0kCFKxJe7Xyd+LNYqNYqoVRP4qG/NOOVF2rxzq/0UKh9+MPPKkE6xVixaRZ//ndDkEtFSXee5RqxiOneM0GTqGhwYcSkxHCvAo2rYt+w/ReZgCu7/QcwwcRsUKximo2EByc3tE/j3JEidB7S4qYHrke878oZNcaH5Cgq9WDVrdF3C0QEwUshU2WRghcFwMUuQAfzn6GDA3kEwtAXyqCSSbD6tlksv2kWLocYr75YA9Ipxvo7p/GVHwxCKqNRp1Pyun2WWE6/G7CUMYsTuQEBBbB2PriVDfMfYFyraIt0aZsd6iIvLINyjPZFulPH/sg8/W/9E2vMTjw6tYyXvwNfGYFczAbcbzQai/2DQgU4hiAIrKstAkEADz89CpmCeWptut+Kx3cOQ+y34XLl/7OwtfaBpJo4k2s8JIGhKwUgCDqkkIbPFOBworkBYgmwbjNj+p//uCTkZ5YhGXqMhZDIKNx8b+xdIL6qMCuLMqssjMWiDgIuFNZUa6CUilGs9+Ivd/fjz/9uCNv+mxlsXIgAkXMFoPc/vZe22SES0+g+VQjrWHSTfqBb6ff/XQH/H+CnB0AsYg1M3dDOKIALx0rgC8q9Oe73/Te0T0OtiZ2Uw3UAkCXaiDMuEBTAAkAqFmG9f6Z87VIX1t4xE9KAUl8o56V+PhVKC5gpPtoyL9ZtngZFETj6hi7qsWz+f1OQ/0+AyKoFEOtJXdvqQlk1CYdNin/960a8tU+P80cLYTxYBAC47Qvxp0NxHQAMppKHGI+gABYIrBUQjYYEgySygUQsgk7FLODNjzAR8lPvFWNuNvIWu3aBkXdJUAFQoULCew/DYIqU0qiDOAkCuPPRSRAiGoNXlPj4D6V45YVakE4xmtc6IvIwgpGIRLwqsWoe4gCCAlggSMUibKiLTE2WS8Ro0Wc//TcabPvrqiYSrRvs8JAinHg71AognUz+Pwg6pJAmm+Y/Syw34KZ7bPjbX3Vh+/8cwJY/HUfzGgeK9W7c8+X4GYqlBTJe+zDwkRK86IqBFjKrqzU4c90WKAwpL1Tg3pV6aHhuGpEsZYVydJuZ7bHNHZO4cqYAn7ypwx3bJiGV0fB5gd/sroHPK0JNqx3Kgnn/n+8agGiUa+Ton4yesqxQ+7DM4MAyQ/Jp1Xqe27CxDUKmnWk0MIyBYAEsIIKtgNXVRXhkQ1XeLH4gtNKuZZ0DVU1O2KckOHO4CBQF/P5/VaH700KoNF489K2+kPdmq4VZMFz763z6/yxVHMcBBAWwwFhdrcF9K8txZ2spxFlq+5UspQUyEGBkIoj5WMDHfyjB/91bjjOHtZApKGz/h+soqwlNo82VBcDKywV87QAEw3U+gKAAFhhSsQhL9AW5FiMqMokIxUHNPNbeMQNtmQfjQ3J88scSiCU0vvo/BlG3NHTxS0QiaHNgySikYhQpufGCZeLQz84XXMcBBAUgwCnBHXfEEuD2hxkrgCBo/OnOYbRuiPSpyzXyrEzBiUZ5jEBgqpQVZuczFKtkEROOMkEIAgpwSlmhDJeDCik33TcFy6AcS9Y6sPaOmajvybR9diZUaORxZ/Clcp5sUVmkgGmCm54PggIQ4JTwllsyBY2O74zGfU9tDhUAV6m7VWmOAEvrd3GoAAQXQIBT9IWpBdZkYlFWouexKC2QZxxMJQgiK1N8WLjcCRAUgACnyCSilAJrVVpl1oaYREMsIqI2Ck2FUrUsq1mMpWoZZBxVHAoKQIBzwltvxyOX5j9Lpgk81Vn+DCIRgYoibqwmQQEIcM5CUwCxUoKTJZvmP9e/U1AAApyTbNGKSibJSQJQOJnEIAhk1/9nqSriRnEKCkCAcyqKFEkl9tQU56aFWThaVfTKwGQoVkuhzOIgFhY9R9OfBAUgwAvLKwsTHpMP5j9LutuBNVnc/guGq7ZjggIQ4IWl5YUJtwNri3Pfx4ClXpeeLFwX52SbjBKBDAaDFsBW/7cbjUZjUhOCBG58ChUS1BQrIsaBs2jkkqy0AE+WplI1jlxNfZpPLvx/LsnUAngMgI5tCBpjhoDAImVZRWw3oLIo98G/YAoVkpSr+YpVspy3YsuUjKQPmxnQBOAX8Y7Pp9mAmZDvMuaLfBqahtvljDoctEDN/wzAVCmWuNHnn8EXbzYgS42Kzuln4OLvzIn6MhgMTQCsRqPRFO+4fJoNmCn5LmO+yLduWoRLo6FFQAQILK1W5Y2MLG0qLS5PDQa+jzb3L5jl9Xro9YmDnXyS6TXkajRYh9Fo/GZGkgjckCyvKIxQAMsqCqDIn0ZGAUoKZCm13MrVDgCXZDwazGAwdBiNxn/0f731RpgZKMAdVVpFYFEppWJsbi1Fi74g78x/lqZSNboGE5vWGoUUBXkUxEyXTHcBtgJ40WAw/J3/JWEXQCCCZRWFmLC7cWdrKVQ5SJpJhWQVQPUC3/5jyTQIeBBAM0eyCNygtNVp865/YSwqi+RQySRwJCi3X1GpyY5APCMkAgnwzkJZ/ABT259ozFqdTrXgE4BYBAUgIBBGc1n86P+mxuIsScI/ggIQEAijRquANEaufWOpmrNGovmAoAAEBMKQiEWoL47MCiRA4OaGG+fpDwgKQEAgKrc1arC2pijktaYydUrNThYCggIQEIiCiCBwR0sp7llRDqlYdEM+/QGhLbiAQFxaywtQopbhisWeF92LuEZQAAICCSgpkOGWAl3iAxcgggsgILCIERSAgMAiRlAAAgKLGEEBCAgsYgQFICCwiBEUgIDAIkZQAAICixhBAQgILGIEBSAgsIgRFICAwCJGUAACAouYjGsB/I1BAeBuYTSYgMDCIiMLwGAwbACz8A8C2OAfECIgILBAyLQr8BkAZ/xDQk2JJgP96Ec/yuTXCQgIcAxX5cAGAHGbqe/evXvhtIYVEFgkEDRNxz0gydFgMBgMvwDwATspWEBAIP/JaDSYwWB4EcA1/zE2ADdm1wQBgRuUhBZAPPxBPzbw96gwIFRAYGGRkQIQEBBY2AiJQAICi5i8agpqMBg6wMQSNrAjx/MJf9LT3xqNxrtzLUs0/NuxbGLWxnxMzFpIiWMGg+HFfJTRYDBMATABOJipfHljAfgXPztx2BZ0o+QN4TsfechjAHTsTkyMHZycsZASx/z3X77K96jRaGzjQjnlkwWwEcBr/q9NADYAyPcFl1eE7dg0AfhFrmSJRqqJY7nCr5jyUjY/WoPB0MTF9csbCwCANuz7kpxIcQPgv4Gt+brAkETiWI7hZHHxiA6A1Z97kxH5pACEPALu6MjnLVm/C6Bl3b58wmAwbM13V89oNO41Go02MK5yRtcwnxTAacxbAU0APsihLAsWg8HQwQZQ8y2OYjAYXgyKS+SrwrcaDIat/oXV5I9b5A0Gg+EbQYt+MtPz5Y0C8AeumtibNh+1sP/CG/LxyQUEFvyLBoOh02AwdOZanij8AoDJL6c2XpZprjAajWf8954OkW5pPvB7BAXJM029FxKBBAQWMXljAQgICGQfQQEICCxiBAUgILCIERSAgMAiRlAAAgKLGEEBCAgsYgQFICCwiBEUgIDAIub/B+UgPARSvJmMAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 288x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Test points every 0.1 between 0 and 5\n",
    "test_x = torch.linspace(0, 5, 51)\n",
    "\n",
    "# Get into evaluation (predictive posterior) mode\n",
    "model.eval()\n",
    "likelihood.eval()\n",
    "\n",
    "# The gpytorch.settings.fast_pred_var flag activates LOVE (for fast variances)\n",
    "# See https://arxiv.org/abs/1803.06058\n",
    "with torch.no_grad(), gpytorch.settings.fast_pred_var():\n",
    "    # Make predictions\n",
    "    observed_pred = likelihood(model(test_x))\n",
    "\n",
    "    # Initialize plot\n",
    "    f, ax = plt.subplots(1, 1, figsize=(4, 3))\n",
    "\n",
    "    # Get upper and lower confidence bounds\n",
    "    lower, upper = observed_pred.confidence_region()\n",
    "    # Plot training data as black stars\n",
    "    ax.plot(train_x.numpy(), train_y.numpy(), 'k*')\n",
    "    # Plot predictive means as blue line\n",
    "    ax.plot(test_x.numpy(), observed_pred.mean.numpy(), 'b')\n",
    "    # Shade between the lower and upper confidence bounds\n",
    "    ax.fill_between(test_x.numpy(), lower.numpy(), upper.numpy(), alpha=0.5)\n",
    "    ax.set_ylim([-3, 3])\n",
    "    ax.legend(['Observed Data', 'Mean', 'Confidence'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
